# -*- coding: utf-8 -*-
"""
Creating a normalised model from the original model

@author: obiri
"""
import mpctools as mpc
import numpy as np
import casadi as cs
import time

##16 states, 21 parameters. 
#nb: some of the remaining parameters are fixed measurements, one is state dependent

lb = 0.1 # lower bound of normalized model: lb*xs
ub = 1.5 # upper bound of normalized model: ub*xs


width = ub-lb

# xs = np.array([3.27367434e+10, 8.28831113e+10, 2.71741538e+03, 5.78604334e+00,
#         6.37656737e+01, 4.83137928e+00, 6.23866477e+01, 3.00039503e+03,
#         2.61960613e+09, 6.63233551e+09, 2.17448574e+03, 4.63001301e+00,
#         5.10255242e+01, 3.86608729e+00, 4.99220220e+01, 2.97694060e+03,
#         3.70349577e+01, 1.50000000e+01, 4.99220218e+01])

xs = np.array([3.27367434e+10, 8.28831113e+10, 2.71741538e+03, 5.78604334e+00,
        6.37656737e+01, 4.83137928e+00, 6.23866477e+01, 
        2.61960613e+09, 6.63233551e+09, 2.17448574e+03, 4.63001301e+00,
        5.10255242e+01, 3.86608729e+00, 4.99220220e+01, 
        3.70349577e+01, 4.99220218e+01])



##THESE PARAMTERS ARE THE NOMINAL PARAMETERS. WE ADD SOME BIAS TO REPRESENT THE TRUE PARAMETERS.
##THEY ARE THEN AUGMENTED AS EXTRA STATES, AND ALL THE AUGMENTED STATES ARE SCALED AND USED FOR ESTIMATION
# p_n = np.array([1.76, 9.6e-03/60, 0.75, 0.075, 28.48, 171.76, 2, 0.45, 2.0, 
#                 2.6e8, 8e8, 4, 0.06/60, 6.59e-10/60, 5e5, 1560.0, 1.244, 4e2])

num_est_param = 18
pn=np.zeros((num_est_param))
pn[0]= 1.76
pn[1]= 9.6e-03/60  #fix  #0.00016
pn[2]= 0.75
pn[3]= 0.075
pn[4]= 28.48
pn[5]= 171.76
pn[6]= 2.0 #fix
pn[7]= 0.45
pn[8]= 2
pn[9]= 2.6e8
pn[10]= 8e8
pn[11]= 4
pn[12]= 0.06/60      #0.001
pn[13]= 6.59e-10/60  #1.0983333333333334e-11
pn[14]= 5e5
pn[15]= 1560.0
pn[16]= 1.244
pn[17]= 4e2

# # ##Generating the true parameters for the plant by adding some bias
# np.random.seed(5)
# q=np.zeros((num_est_param))
# # q = np.random.normal(loc=q, scale=std_dev)    
# q[0] = np.random.normal(loc=pn[0], scale=0.0001)   
# q[1] = np.random.normal(loc=pn[1], scale=0.000001)   
# q[2] = np.random.normal(loc=pn[2], scale=0.00001) 
# q[3] = np.random.normal(loc=pn[3], scale=1e-6)   
# q[4] = np.random.normal(loc=pn[4], scale=1e-3) 
# q[5] = np.random.normal(loc=pn[5], scale=1e-3)   
# q[6] = np.random.normal(loc=pn[6], scale=1e-5) 
# q[7] = np.random.normal(loc=pn[7], scale=1e-5) 
# q[8] = np.random.normal(loc=pn[8], scale=1e-4) 
# q[9] = np.random.normal(loc=pn[9], scale=1e4) 
# q[10] = np.random.normal(loc=pn[10], scale=1e4) 
# q[11] = np.random.normal(loc=pn[11], scale=1e-3) 
# q[12] = np.random.normal(loc=pn[12], scale=1e-8) 
# q[13] = np.random.normal(loc=pn[13], scale=1e-11) 
# q[14] = np.random.normal(loc=pn[14], scale=1e2) 
# q[15] = np.random.normal(loc=pn[15], scale=1e-1) 
# q[16] = np.random.normal(loc=pn[16], scale=1e-3) 
# q[17] = np.random.normal(loc=pn[17], scale=1e-3) 

factor = 1.01*np.ones(18)
factor[0] = 1.1
factor[3] = 1.1
factor[12] = 1.1
factor[13] = 1.1
factor[15] = 1.1


p_n =factor*pn

# xs_aug = np.concatenate((xs,p_n))

def reactor_tank(x_aug, u, w):   
    
    """
    An integrated model including upstream and buffer tank only.
    :param x: State vector (19,1)
    :param u: Input vector (8,1)
    :return: ODEs (19,1)
    """
    dxdt = cs.SX.zeros(x_aug.shape[0])
    
    x = x_aug[0:16]
    p = x_aug[16:]
       
    #new state indexing
    # Xv1 = x[0]   # concentration of viable cells in reactor(cell/L)
    # Xt1 = x[1]   # total concentration of cells in reactor(cell/L)
    # GLC1 = x[2]  # glucose concentration in reactor(mM)
    # GLN1 = x[3]  # glutamine concentration in reactor(mM)
    # LAC1 = x[4]  # lactate concentration in reactor(mM)
    # AMM1 = x[5]  # ammonia concentration in reactor(mM)
    # mAb1 = x[6]  # mAb concentration in reactor(mg/L)
    # Xv2 = x[7]    # concentration of viable cells in separator(cell/L)
    # Xt2 = x[8]    # total concentration of cells in separator(cell/L)
    # GLC2 = x[9]  # glucose concentration in separator(mM)
    # GLN2 = x[10]  # glutamine concentration in separator(mM)
    # LAC2 = x[11]  # lactate concentration in separator(mM)
    # AMM2 = x[12]  # ammonia concentration in separator(mM)
    # mAb2 = x[13]  # mAb concentration in separator(mg/L)
    # # temperature as a state -- Ben
    # T = x[14]     #temperature of bioreactor mixture  (d C)
    # # buffer tank
    # c = x[15]  # Concentration of mAb in the tank (mg/L)


    #defining new variables for the normalised states
    x1 = x[0]*(width*xs[0])+lb*xs[0];      x2 = x[1]*(width*xs[1])+lb*xs[1];      x3 = x[2]*(width*xs[2])+lb*xs[2]; 
    x4 = x[3]*(width*xs[3])+lb*xs[3];      x5 = x[4]*(width*xs[4])+lb*xs[4];      x6 = x[5]*(width*xs[5])+lb*xs[5]; 
    x7 = x[6]*(width*xs[6])+lb*xs[6];      x8 = x[7]*(width*xs[7])+lb*xs[7];      x9 = x[8]*(width*xs[8])+lb*xs[8];
    x10 = x[9]*(width*xs[9])+lb*xs[9];     x11= x[10]*(width*xs[10])+lb*xs[10];   x12 = x[11]*(width*xs[11])+lb*xs[11];
    x13 = x[12]*(width*xs[12])+lb*xs[12];  x14= x[13]*(width*xs[13])+lb*xs[13];   x15 = x[14]*(width*xs[14])+lb*xs[14];
    x16 = x[15]*(width*xs[15])+lb*xs[15];  
    
    p1 =  p[0]*(width*p_n[0])+lb*p_n[0];   p2 =  p[1]*(width*p_n[1])+lb*p_n[1];     p3 =  p[2]*(width*p_n[2])+lb*p_n[2];
    p4 =  p[3]*(width*p_n[3])+lb*p_n[3];   p5 =  p[4]*(width*p_n[4])+lb*p_n[4];     p6 =  p[5]*(width*p_n[5])+lb*p_n[5];    
    p7 =  p[6]*(width*p_n[6])+lb*p_n[6];   p8 =  p[7]*(width*p_n[7])+lb*p_n[7];     p9 =  p[8]*(width*p_n[8])+lb*p_n[8]; 
    p10 = p[9]*(width*p_n[9])+lb*p_n[9];   p11 = p[10]*(width*p_n[10])+lb*p_n[10];  p12 = p[11]*(width*p_n[11])+lb*p_n[11];
    p13 = p[12]*(width*p_n[12])+lb*p_n[12];p14 = p[13]*(width*p_n[13])+lb*p_n[13];  p15 = p[14]*(width*p_n[14])+lb*p_n[14];
    p16 = p[15]*(width*p_n[15])+lb*p_n[15];p17 = p[16]*(width*p_n[16])+lb*p_n[16];  p18 = p[17]*(width*p_n[17])+lb*p_n[17];

    
    #inputs remained the same as defined above
    # x8 = xs[7]; 
    # x16 = xs[15];
    # x18 = xs[17];
    
    # Inputs of upstream model:
    F_in = u[0]    # inlet flow rate (L/min)
    F_1 = u[1]     # outlet flow rate of bioreactor (L/min)
    F_r = u[2]     # u[1]-u[3] # u[2]      # recycle flow rate (L/min)
    F_2 = u[3]     # Outlet flow rate of separator (L/min)
    GLC_in = u[4]  # inlet glucose concentration (mM)
    GLN_in = u[5]  # inlet glutamine concentration (mM)
    Tc = u[6]      # coolant temperature (d C)
    # Inputs of buffer tank:
    F_out_bf = u[7]  # Outlet flow rate of the buffer tank (L/min)  
    F_in_bf = u[3]   # Inlet flow rate of the buffer tank (L/min)
    #c_in_bf = x[14] # Inlet concentration of mAb (mg/L)
    c_in_bf = x14    ##changed to the scaled value of x[14]
    

    #18 selected Parameters as augmented states:
    # K_damm = p[0]#1.76   # ammonia constant for cell death (mM)
    # K_dgln = p[1]#9.6e-03/60   # constant for glutamine degradation (min^(-1))
    # K_glc =  p[2]#0.75   # Monod constant for glucose (mM)
    # K_gln =  p[3]#0.075  # Monod constant for glutamine (mM)
    # KI_amm = p[4]#28.48  # Monod constant for ammonia (mM)
    # KI_lac = p[5]#171.76   # Monod constant for lactate (mM)
    m_glc  = 4.9e-14/60   # maintenance coefficients of glucose (mmol/cell/min)
    # n      =   p[6]#2    # constant represents the increase of specific death as ammonia concentration increases (-)
    # Y_ammgln = p[7]#0.45 # yield of ammonia from glutamine (mmol/mmol)
    # Y_lacglc = p[8]#2.0  # yield of lactate from glucose (mmol/mmol)
    # Y_Xglc = p[9]#2.6e8  # yield of cells on glucose (cell/mmol)
    # Y_Xgln = p[10]#  8e8   # yield of cells on glutamine (cell/mmol)
    alpha1 = 3.4e-13/60  # constants of glutamine maintenance coefficient (mM L/cell/min)
    # alpha2 = p[11]# 4  # constants of glutamine maintenance coefficient (mM)
    mu_max = (0.0016*x15 - 0.0308)*5/60 # 0.058 * 5  ##x[16] changed to the scaled value of x16
    # mu_dmax = p[12]# 0.06/60
    # m_mabx =  p[13]# 6.59e-10/60  # constant production of mAb by viable cells (mg/Cell/min)

    # Add new parameters -- Ben
    # Delta_H = p[14]# 5e5   # 4.4e4# 1e-8
    # rho =     p[15]#1560.0      # density of mixture is assumed to be density of glucose [g/L]
    # cp  =     p[16]#1.244      # specific heat capacity of mixture is assumed to be that of glucose [J/g/dC]
    # U   =     p[17]#4e2
    T_in = 37.0   ##not included in the augmented system because temperature can easily be measured
    
    ####---these parameters were excluded because they are very easy to know, they willnot change even after a very long time----####
    # Parameters of buffer tank
    D = 5  # Diameter of the tank (dm)
    L = D  # Diameter of the tank (dm)
    Ac= np.pi*(D/2)**2  # Cross-sectional area (dm^2)
    V = Ac*L  # Volume of the tank (L)
    
    
    # Assumptions: 92% cell recycle rate and 20% mab retention
    Xvr = 0.92 * x1 * F_1 / F_r  # concentration of viable cells in recycle stream(cell/L)
    Xtr = 0.92 * x2 * F_1 / F_r  # total concentration of cells in recycle stream(cell/L)
    GLCr = 0.2 * x3 * F_1 / F_r  # glucose concentration in recycle stream(mM)
    GLNr = 0.2 * x4 * F_1 / F_r  # glutamine concentration in recycle stream(mM)
    LACr = 0.2 * x5 * F_1 / F_r  # lactate concentration in recycle stream(mM)
    AMMr = 0.2 * x6 * F_1 / F_r  # ammonia concentration in recycle stream(mM)
    mAbr = 0.2 * x7 * F_1 / F_r  # mAb concentration in recycle stream (mg/L)
    
    
    # ODEs
    f_lim = (x3 / (p3 + x3)) * (x4 / (p4 + x4))
    f_inh = (p6 / (p6 + x5)) * (p5 / (p5 + x6))
    mu = mu_max * f_lim * f_inh
    mu_d = p13 / (1 + (p1 / x6) ** p7)
    
    dxdt[0]= (mu * x1 - mu_d * x1 - (F_in / 3.00039503e+03) * x1 + (F_r / 3.00039503e+03) * (Xvr - x1) +w[0]) /(width*xs[0])
    dxdt[1] = (mu * x1 + (F_r / 3.00039503e+03) * (Xtr - x2) - (F_in / 3.00039503e+03) * x2 +w[1]) /(width*xs[1])
    
    Q_glc = mu / p10 + m_glc
    dxdt[2] = (-Q_glc * x1 + (F_in / 3.00039503e+03) * (GLC_in - x3) + (F_r / 3.00039503e+03) * (GLCr - x3) +w[2]) /(width*xs[2])
    
    m_gln = (alpha1 * x4) / (p12 + x4)
    Q_gln = mu / p11 + m_gln
    dxdt[3] = (-Q_gln * x1 - p2 * x4 + (F_in / 3.00039503e+03) * (GLN_in - x4) + (F_r / 3.00039503e+03) * (GLNr - x4) +w[3]) /(width*xs[3])
    
    Q_lac = p9 * Q_glc
    dxdt[4] = (Q_lac * x1 - (F_in / 3.00039503e+03) * x5 + (F_r / 3.00039503e+03) * (LACr - x5) +w[4]) /(width*xs[4])
    
    Q_amm = p8 * Q_gln
    dxdt[5] = (Q_amm * x1 + p2 * x4 - (F_in / 3.00039503e+03) * x6 + (F_r / 3.00039503e+03) * (AMMr - x6) +w[5]) /(width*xs[5]) 
    dxdt[6] = (p14 * x1 - (F_in / 3.00039503e+03) * x7 + (F_r / 3.00039503e+03) * (mAbr - x7) +w[6]) /(width*xs[6])
    
    # dxdt[7] = (F_in + F_r - F_1)/(width*xs[7])
    # dxdt[7] = 0
    
    dxdt[7] = ((F_1 / 2.97694060e+03) * (x1 - x8) - (F_r / 2.97694060e+03) * (Xvr - x8) +w[7]) /(width*xs[7])
    dxdt[8] = ((F_1 / 2.97694060e+03) * (x2 - x9) - (F_r / 2.97694060e+03) * (Xtr - x9) +w[8]) /(width*xs[8])
    dxdt[9] = ((F_1 / 2.97694060e+03) * (x3 - x10) - (F_r / 2.97694060e+03) * (GLCr - x10) +w[9]) /(width*xs[9])
    dxdt[10] = ((F_1 / 2.97694060e+03) * (x4 - x11) - (F_r / 2.97694060e+03) * (GLNr - x11) +w[10]) /(width*xs[10])
    dxdt[11] = ((F_1 / 2.97694060e+03) * (x5 - x12) - (F_r / 2.97694060e+03) * (LACr - x12) +w[11]) /(width*xs[11])
    dxdt[12] = ((F_1 / 2.97694060e+03) * (x6 - x13) - (F_r / 2.97694060e+03) * (AMMr - x13) +w[12]) /(width*xs[12])
    dxdt[13] = ((F_1 / 2.97694060e+03) * (x7 - x14) - (F_r / 2.97694060e+03) * (mAbr - x14) +w[13]) /(width*xs[13])
    # dxdt[15] = (F_1 - F_2 - F_r) /(width*xs[15])
    # dxdt[15] = 0
    
    # Temperature dynamics -- Ben
    # Average mass of an mAb cell is 150 kDa = 2.4908084e-19 g
    dxdt[14] = ((F_in / 3.00039503e+03) * (T_in - x15) + (p15 / (p16 * p17)) * mu * x1 * 2.4908084e-19 + (p18 / (3.00039503e+03 * p16 * p17)) * (Tc - x15) +w[14]) /(width*xs[14]) 
    
    
    # ODEs of buffer tank:
    dhdt = 1/Ac*(F_in_bf-F_out_bf)
    dcdt = F_in_bf/(Ac*1.50000000e+01)*(c_in_bf-x16)
    # dxdt[17] = dhdt /(width*xs[17])
    # dxdt[17] = 0
    dxdt[15] = (dcdt +w[15]) /(width*xs[15])
                                                                                                                                                         
    #dxdt = [dxdt1, dxdt2, dxdt3, dxdt4, dxdt5, dxdt6, dxdt7, dxdt8, dxdt9, dxdt10, dxdt11, dxdt12, dxdt13, dxdt14, dxdt15, dxdt16, dxdt17, dxdt18]
    return dxdt #(returns scaled values since this is now the scaled version of the ode )



Ny = 7
Nx = 16
Nx_aug = 16 + 18

def measurement_xy_model(x_aug):            
    x = x_aug[0:16]
    
    C_matrix = np.zeros((Ny, Nx))
    C_matrix[0,7]=1
    C_matrix[1,8]=1
    C_matrix[2,9]=1
    C_matrix[3,10]=1
    C_matrix[4,11]=1
    C_matrix[5,12]=1
    C_matrix[6,15]=1   
    y = np.dot(C_matrix, x)
    return y








